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Abstract 

The explicit analytic solution of the Thomas- Fermi equation thorough 
a new kind of analytic technique, namely the homotopy analysis method, 
was employed by Liao [?] (Appl. Math. Comp. 144, (2003)). However, 
the base functions and the auxiliary linear differential operator chosen were 
such that the convergence to the exact solution was fairly slow. New base 
functions and auxiliary linear operator to form a better homotopy are the 
main concern of the present paper. It is known that proper choice of base 
functions and auxiliary operator is extremely significant in gaining the 
exact solution in order to reduce the computational cost. The proposed 
homotopy here not only greatly reduces the computational efforts by at 
least doubling the convergence of the homotopy series, but also enlarges the 
convergence region of the homotopy series as compared with that of Liao 
[?]. Fade approximants to the obtained solutions increase the accuracy 
even to a higher degree. To support this, the explicit analytical expressions 
obtained using the proposed approach are compared with the numerically 
computed ones and those of Liao [?] . 
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1. Introduction 

Nonlinear phenomena are encountered in all areas of sciences and engineering. 
Even though the study of nonlinear equations is of great importance to many sci- 
entific researchers in various fields, it is very difficult to solve nonlinear problems 
and, in general, it is often more costly to get an analytic approximation than a 
numerical one to a given nonlinear problem [?]. 
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Since the Thomas-Fermi equation is one typical not easy-to-find exact solu- 
tion, the explorations of the solution have taken considerable attraction of many 
authors. The analytic approximations of the Thomas-Fermi equation were pro- 
posed by some different techniques such as the variational approach [?], the 5- 
expansion method [?, ?], the decomposition method [?, ?], see also [?, ?, ?]. Since 
the techniques used in these papers were mostly perturbative and limited, the 
solutions are obliged to be valid for certain regions with a decreasing accuracy. 

In an aim to search for new and more effective analytical tools, Liao in 1992 
[?, ?] proposed the homotopy analysis method which deforms a difficult nonlin- 
ear problem into easy linear counterparts. This is achieved by introducing an 
auxiliary parameter in the construction of an homotopy, which can provide a 
convenient way to control the convergence of the approximation series and adjust 
the convergence regions if necessary. A series of nonlinear problems were later 
attacked with the use of homotopy analysis method [?, ?, ?, ?, ?, ?] and [?, ?]. 
An elegant, simple and explicit analytic solution to the Thomas- Fermi equation 
was presented in [?]. One main drawback of the homotopy analysis method is 
that it's convergence to the exact solution becomes too time-consuming as the 
parameters involved in the equations are not treated rationally. This seems to 
be the reason why Liao in [?] was unable to obtain higher-order approximations 
for the solution of the Thomas-Fermi equation, but presented only the homotopy 
approximations with a percentage accuracy of at most first order. On the other 
hand, this order is inadequate to demonstrate convergence of the method to the 
true solution. 

We in the present paper use homotopy analysis technique of Liao for the ana- 
lytic calculation of Thomas- Fermi equation. Motivated by the study [?], we aim 
to show that the computational task during the implementation of the homotopy 
analysis method can be further reduced if the Thomas-Fermi equation is treated 
in a more rational manner. This is accomplished here by modifying the original 
approach of Liao [?]. Knowing that the choice of base functions and auxiliary 
linear differential operator is very important in the homotopy technique, see for 
example [?], new base functions and auxiliary linear differential operator are pre- 
sented here. Unlike those of [?] , using these new parameters it is shown here that 
the convergence of the homotopy can be greatly accelerated, even greater by the 
homotopy Fade approximants. In addition to this, the region of convergence is 
much improved with the present approach. 

The following strategy is adopted in the rest of the paper. In §5] the basis 
of homotopy analysis method is laid out with an application to the nonlinear 
Thomas-Fermi equation. Analytic expressions for the solution are derived and 
compared with the numerical and previously published results in ^2 Finally 
conclusions follow in ^13 
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2. The Homotopy Analysis Method 



The Homotopy analysis method was first proposed by the Chinese mathemati- 
cian Liao [?]. This method is based on the homotopy and has several advantages. 
To underline, firstly its validity does not depend upon whether or not nonlin- 
ear equations under consideration contain small or large parameters, hence it 
can solve more of strongly nonlinear equations than the perturbation techniques. 
Secondly, it provides us with a great freedom to select proper auxiliary linear 
operators and initial guesses so that uniformly valid approximations can be ob- 
tained. Thirdly, it gives a family of approximations which are convergent in a 
larger region. Liao successfully applied the homotopy analysis method to solve 
some nonlinear problems in mechanics. For example, Liao in [?] gave a purely 
analytic solution of 2D Blasius's viscous flow over a semi-infinite flat plate, which 
is uniformly valid in the whole physical region. 

The essential idea of this method is to introduce a homotopy parameter, say 
p, which varies from to 1. At p = 0, the system of equations usually has been 
reduced to a simplified form which normally admits a rather simple solution. As p 
gradually increases continuously toward 1, the system goes through a sequence of 
deformations, and the solution at each stage is close to that at the previous stage 
of the deformation. Eventually at p = 1, the system takes the original form of the 
equation and the final stage of the deformation gives the desired solution. Here 
we apply it to solve the nonlinear Thomas-Fermi problem. Consider a differential 
equation used to calculate the electrostatic potential in the Thomas-Fermi atom 
model [?, ?], called the Thomas- Fermi equation with the boundary conditions 



which describes the spherically symmetric charge distribution about a many elec- 
tron atom. Equation (12. ip is not amenable to exact treatment and, therefore, 
approximate techniques must be resorted to. There exists neither linear terms 
nor small or large parameters in equation (12.11) . so the standard perturbation 
methods cannot be applied directly. Due to the fact that the homotopy analy- 
sis method requires neither a small parameter nor a linear term in a differential 
equation, one possibility to approximately solve equation (12. ip is by means of the 
homotopy analysis method. 

We first rewrite the original equation (12. ip in the form 



The essence to approximate a problem is to represent its solution by means 
of a complete set of base functions. Considering the boundary conditions in (12. 2p 
and the physical meaning of u, Liao [?] chosen the set of base functions 






(2.2) 




(2.3) 



3 



As Liao [?] discusses in detail on the free falling problem, the selection of 
different sets of functions is possible all of which generate the same solution. 
Besides, the homotopy analysis method provides us with freedom to choose the 
initial guess and the auxiliary linear operator so that one can represent the so- 
lution of the Thomas-Fermi equation by distinct set of base functions. However, 
some set of base functions might cause a very slow convergence to the true so- 
lution. Therefore, if a proper choice is made while selecting a set of function 
a better convergence can be achieved which may save both time and resources 
while computing. Thus, we here take a different set of base functions from (12.31) 
in an aim to speed up the convergence to the exact solution and instead of (12.31) 
choose 



{a{a + (3x)-^"'\m>l, a,P,jeR} (2.4) 
to represent the solution u{x) of (12.21) in the form 

°° a 

u(x) = c„- — — , (2.5) 

where c^'s are coefficients. Notice that if we set a = /3 = 7 = 1 in (I2.4ti2.5l) . the 
base functions of Liao ([?]) in equation (12. 3p are obtained. Therefore, our base 
functions are somewhat more general. The key to choose the parameters in (12.41) 
is that a better convergence rate can be captures as compared to the case of Liao 
[?]. It will be clear soon that some certain selections will indeed lead to high 
savings in the computational time of the approximate solution of Thomas- Fermi 
equation by accelerating its convergence. 

Equation (12.41) provides us with the rule of solution expression. This rule 
is important in the frame of the homotopy analysis method. Considering the 
initial conditions in (12.21) and the rule of solution expression above described, it 
is obvious that 

a 

is a good initial guess for u{x) satisfying the boundary conditions in (12. 2p exactly. 
We next choose here 

^ Px + a d'^ ^ ^ d 
a + J dx'^ dx 

as our auxiliary linear differential operator, which was found to be quite efficient 
for the consideration of the present nonlinear problem, which has the property 
that L{Ci{a + (3x)~"' + C2) = 0, with Ci and C2 arbitrary constants. We are 
now at the stage of constructing the zeroth order deformation equation system 
associated with the Thomas- Fermi problem given by (12.21) . which is also called a 
family of differential equations (viewing p as a parameter) 

{l-p)L{u) = phN{u), u{0,p)-l = u{oo,p) = 0. (2.7) 
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It should be noted that the nonhnear differential operator in equation fl2.5p is 
given as in (12.21) . Moreover, the parameter h is an auxiliary non-zero parameter 
to adjust the convergence rate of the perturbation series, which was found to be 
— 1/2 < /i < in [?] for the convergence. However, in our case we will show that 
the region of h for the convergence can be extended. Obviously when p = and 
p = 1, we have respectively 

0) = Uo(x), 1) = m(x). (2.8) 

Hence the process of giving an increment to p from to 1 is the process of 
u{x,p) varying continuously from the initial guess Uo{x) to the final solution u{x). 
This kind of continuous variation is called deformation in topology so that we 
call system (12. 7p the zeroth order deformation equation. Next, differentiating 
(12. 7p successively and eventually imposing at p = 0, the kth-order deformation 
equations follow as 

L{uk) = L{uk-i)+hRk, Mfc(O) =Mfc(cx)) =0. (2.9) 

The function i?^ on the right-hand side of (12.91) is given by 

1 d^-'^N 

- (A;-l)! 9j9^-l 'P=°- 

Finally, a straightforward Taylor expansion of u{x,p) at the point p = and 
eventually imposing the series at p = 1 gives the solution of system (12.21) in the 
form 

oo 

u{x) = uo{x) + ^Uk{x), (2.10) 

k=l 

for which we presume that the initial guesses Uq to u, the auxiliary linear oper- 
ator L and the non-zero auxiliary parameter h are all so properly selected that 
the deformations u{x,p) are smooth enough and their kth-order derivatives with 
respect to p in equation (12.91) exist and are given by Uk = ■^|^|p=o- It is clear 
that the convergence of Taylor series at p = 1 is a prior assumption here so that 
the system in (12.101) holds true. The formulae in (12.101) provide us with a di- 
rect relationship between the initial guesses and the exact solutions. Moreover, a 
special emphasize should be placed here that the kth-order deformation system 
(12. 9p is a linear differential equation system with the auxiliary linear operator L 
whose fundamental solution is known as aforementioned. Eventually, we obtain 
the result at the Nth-order approximation as 

N 

uix) = (2.11) 

k=0 
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3. Results and discussion 



In this section analytic approximate solutions corresponding to system f l2.2p is 
presented. Solutions obtained from the homotopy analysis method are compared 
with those obtained from the full numerical computations and the available so- 
lution of [?]. 

Note that the analytic solution presented in (12. lip contains the auxiliary 
parameter /i, which can be employed to control the convergence of approximations 
and adjust convergence regions when necessary. As claimed before, the proper 
choice of parameter a, /5, 7 can yield an extended convergence region as compared 
to the case of Liao [?] . 




-1.0 -0,8 -0.6 -0.4 -0.2 0.0 -1.0 -0.8 -0.6 -04 -0.2 0.0 

h h 
(a) (b) 

Figure 1: The h— curves of ti'(O) obtained from the 20th-order homotopy ap- 
proximation solution of the problem (12. 2p . (a) Case of Liao [?] (a = /3 = 7 = 1), 
(b) the present case (a = 3/4, /5 = 7 = 1). 



To illustrate this, the influence of h on the convergence of the solution series 
are given in figures [It^a-b). The /i— curves are drawn for 'u'(O) obtained from the 
20th-order homotopy analysis approximation. Figure [TJ^a) depicts the case of [?] 
(not given in [?]) and[T]^b) depicts our case. As found by Liao in [?], when h is 
restricted to —1/2 < h < 0, the series in (12. lip converges in the whole region 
< X < 00, which is indeed the case as shown in figure [T](a). On the other hand, 
for the parameters a = 3/4, /3 = 7 = 1, we are able to extend the convergence 
region of h to —9/10 < /i < 0, as shown in figure [Tt^b). This is important since 
in a larger region of convergence a better convergence rate can be achieved for 
the series. In fact, it is easy to see that in order to have a good approximation 
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h has to be chosen around —4/5. When h = —4/5, the analytic result at the 
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Figure 2: Comparison of the analytic result of the Thomas- Fermi equation with 
the numerical result. Solid line denotes analytic result at the 40th-order approx- 
imation for h = —4/5 and dashes denote the numerical result. 

40th-order of approximation agrees well with the numerical result, as shown in 
figure [21 As compared with the 60th-order of approximation of Liao [?] (see fig- 
ure 2 of Liao [?]) our analytic solution at a lower order is clearly seen to better 
approach the exact solution. To explain the reason deeply, and in fact to disclose 
the difference between our solution and that of Liao [?], we list in table [1] the 
values of initial slope u'{0) (which is used to obtain the energy of a neutral atom 
in the ThomasFermi model, see [?]) and curvature u"{0) respectively. The error 
between the exact and approximate m'(0) is also given in table [H Obviously, the 
error decreases as the order of approximation increases, in both our case and the 
case of Liao [?] . However, more importantly, the convergence of our series (12.111) 
is seen to be more than twice as the convergence of the iterations of Liao [?]. This 
justifies the claim that our base functions indeed approximate the exact solution 
of Thomas- Fermi equation at a better convergence rate. This as a result will save 
us both in terms of computational time and resources. Moreover, owing to equa- 
tion (12.21) . it holds u"{0) oo as x 0^. The approximations of u"{0) obtained 
from the homotopy analysis method when h = —3/4 are listed in tabled] show 
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that ^"(O) of the analytic solution (12.21) indeed tends to infinity. Again, the trend 
is observed to be twice as high as the trend of Liao [?]. The remarkable accuracy 
of the results tabulated clearly illustrate the effectiveness and efficiency of the 
proposed approach. Table further reveals as claimed before that the presented 
approximate solutions take less computational time as compared with that of [?]. 





u'iO) 


Error (%) 


m'(0)^ 


u"iO) 


u"(0)^ 


10 


-1.54628 


2.63 


-1.50014 


25.4567 


13.0003 


20 


-1.56597 


1.39 


-1.54093 


46.8426 


23.0819 


30 


-1.57305 


0.94 


-1.55595 


68.1948 


33.1119 


40 


-1.57669 


0.71 


-1.56373 


89.5378 


43.1275 


50 


-1.57891 


0.57 


-1.56848 


110.877 


53.1370 


60 


-1.58040 


0.48 


-1.57168 


132.214 


63.1434 


70 


-1.58171 


0.40 


-1.57399 


151.216 


73.1480 


80 


-1.58303 


0.31 


-1.57572 


173.012 


83.1514 


90 


-1.58424 


0.24 


-1.57708 


196.871 


93.1542 


100 


-1.58515 


0.18 


-1.57816 


224.112 


103.1560 



Table 1: Analytic approximations of u'{0) and u"{0) when h = —3/4, the per- 
centage error occurred for u'{0) and comparison of the results with those of Liao 
[?], denoted by a superscript L. 

As implemented in [?], we can further employ the [m,m] diagonal homotopy 
Pade approximants [?] to the power series of u'{0) in order to gain more accurate 
approximations of the initial slope, as shown in table [2l Note that the error 
decreases with the increase of the degree of the Pade approximants. Comparisons 
with the Pade approximations of Liao [?] once more shows the better accuracy 
obtained from the present approach. 

Pade approximants m'(0) Error (%) u'{0)^ 



[10, 10] 


-1.58030 


0.48933 


-1.51508 


[20, 20] 


-1.58571 


0.14867 


-1.58281 


[30, 30] 


-1.58694 


0.07122 


-1.58606 


[40, 40] 


-1.58752 


0.03469 


-1.58668 


[50, 50] 


-1.58801 


0.00384 


-1.58712 



Table 2: Approximations of the initial slope u'{0) given by the diagonal Pade 
approximants, when h = —3/4 and the percentage error occurred. Comparisons 
with u'{0) obtained by Liao [?] is denoted by the superscript L. 
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4. Concluding remarks 



In this paper, the homotopy analysis method has been apphed to obtain ap- 
proximate analytical solution for nonlinear phenomena governed by the Thomas- 
Fermi equation. The results presented have readily revealed that the approach 
adopted is very effective and convenient. Comparisons of the approximations 
with the published ones of Liao [?] have proven the accuracy and efficiency of 
the proposed approach which rapidly doubles the convergence rate of the homo- 
topy series solution to the exact solution and hence significantly reduces the time 
consumption while evaluating the approximate analytic solutions thorough the 
homotopy analysis method. 

The key to get fast convergence and better accuracy as compared with that 
of Liao [?] has been to select more appropriate different set of base functions and 
auxiliary linear differential operator. Such a rational choice has been shown to en- 
large the region of convergence of the homotopy series and thus yield acceleration 
of the convergence. Diagonal homotopy Fade approximants to the power series 
obtained have been shown to improve some degree the accuracy and convergence 
of the homotopy series to the exact ones. 

The purely explicit analytical solutions obtained here also provide a good 
scientific base for the validation of the numerically computed values using different 
schemes in the literature. In addition to this, the developed approach could be 
expected to be applicable to attain the solutions of highly nonlinear systems 
arising from the applications in engineering and science. 
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